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Abstract 

We present a new numerical method for calculating an evolving 2-D Hele-Shaw in- 
terface when surface tension effects are neglected. In the case where the flow is directed 
from the less viscous fluid into the more viscous fluid, the motion of the interface is 
ill-posed; small deviations in the initial condition will produce significant changes in 
the ensuing motion. This situation is disastrous for numerical computation, as small 
round-off errors can quickly lead to large inaccuracies in the computed solution. Our 
method of computation is most easily formulated using a conformal map from the fluid 
domain into a unit disk. The method relies on analytically continuing the initial data 
and equations of motion into the region exterior to the disk, where the evolution prob- 
lem becomes well-posed. The equations are then numerically solved in the extended 
domain. The presence of singularities in the conformal map outside of the disk intro- 
duces specific structures along the fluid interface. Our method can explicitly track the 
location of isolated pole and branch point singularities, allowing us to draw connec- 
tions between the development of interfacial patterns and the motion of singularities 
as they approach the unit disk. In particular, we are able to relate physical features 
such as finger shape, side- branch formation, and competition between fingers to the 
nature and location of the singularities. The usefulness of this method in studying the 
formation of topological singularities (self-intersections of the interface) is also pointed 
out. 

* Research was supported by the National Aeronautics and Space Administration under NASA 
Contract No. NASl-19480 while the author was in residence at the Institute for Computer Applications 
in Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23681-0001. 
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1 Introduction 


The displacement of a viscous fluid by a less viscous fluid in a Hele-Shavv cell has been 
the subject of intense investigation over the last decade, mainly due to newly discov- 
ered mathematical analogies with dendritic crystal growth, directional solidification, 
and electro-chemical growth. The original motivation behind the pioneering work of 
Saffman k Taylor [36] was the analogy to displacement in porous medium. Reviews 
by Saffman [38], Bensimon, Kadanoff, Liang, Shraiman, k Tang [5], and Homsy [16] 
summarize the state of affairs as of the mid-eighties, while, some of the the more recent 
developments are reviewed by Pelce [31], Kessler, Kopiik, k Levine [22], Howison [21], 
and Tanveer [42] from a range of different perspectives. 

In this paper, we shall limit our investigations to channel flow, although our numer- 
ical method is quite general, and can be applied to other geometries. The advantage 
of this restriction is that our description can be specific. Moreover, the phenomena 
displayed in channel flow have been studied extensively, and are representative of other 
Hele-Shaw flows. 

A steadily advancing flat interface in a channel is unstable to perturbations when 
driven by the less viscous fluid. The perturbations grow into fingers, but the subse- 
quent behavior depends on the relative strength of capillary effects, measured by the 
dimensionless number B - ab 2 /(12/iVa 2 ). Here a is the surface tension coefficient, 6 
is the gap width of the cell, a is the channel width, V is the speed of displacement well 
in front of the interface, and is the viscosity of the more viscous fluid (the viscosity of 
the less viscous fluid is assumed negligible). Numerical computations by Tryggvasan 
k Aref [45, 46], DeGregoria and Schwartz [9], Bensimon [6], and Meiburg and Homsy 
[27] show competition between fingers resulting in the emergence of a single steady 
finger, provided B is greater than about 0.0004 1 (but not greater than B c = 0.025 
otherwise the interface is stabilized by capillary effects). For still smaller B , DeGre- 
goria and Schwartz [9, 10] and Bensimon [6] find that the finger spontaneously splits; 
this can be induced for higher values of B by introducing perturbations at the tip, 

*There is a range in B , depending on the level of noise in an experiment or numerical calculation, 
in which a transition occurs between steadily moving fingers and continual unsteady motion. A typical 
range is 0.0005 > B > 0.0002. We have taken 0.0004 as a representative value. 
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even with small amplitude. Experiments also reveal the emergence of tip splitting and 
side-branching instabilities (Park k Hornsy [28], Tabeling, Zocchi, k Libchaber [41]). 
Bensimon [6] provides numerical evidence supporting heuristic arguments that the size 
of the perturbation triggering instability decreases with B. The computations of Dai 
and Shelley [8] (in the circular geometry) also show great sensitivity to the level of 
numerical precision when the surface tension coefficient is small. For small enough B , 
even noise during experiments can be large enough to set off a pattern of continual 
tip-splitting and finger competition ( Maxworthy [26], Arneodo, Couder, Grasseau, 
Hakim k Rabaud [2]). Indeed, when capillary effects are very small, it appears that 
the pattern is fractal (Maxworthy [26], Kopf-Sill k Homsy [23], Arneodo et. al. [2]). 

Detailed understanding of this unsteady behavior is limited. In particular, the 
numerical work has not produced a clear understanding of the asymptotic trends as 
5 — 0. Unfortunately, the inclusion of surface tension in the usual mathematical 
model makes theoretical studies very difficult. Instead a large body of knowledge has 
been developed for the initial-value problem when 5 = 0. For example, Gustaffson 
[13, 14] has rigorously proved the existence of a solution for finite time starting with 
analytic initial data. Earlier, Galin [11] and Polubarinova-Kochina [32] considered 
the mathematically identical problem of the Darcy model for ground water flow and 
devised analytical techniques to obtain exact solutions for a class of initial conditions. 
These were apparently well-known in the Russian literature (see Hohlov [15] and How- 
ison [21]). Exact solutions due to Saffman [37], Howison [18, 19, 20] and Shraiman 
k Bensimon [39] can be seen as applications of these techniques though these results 
were obtained without know-ledge of the earlier Russian work. Howison [21] summa- 
rizes the relation between the different techniques. Within the class of known exact 
solutions, there are finger patterns that exist for all times and exhibit behavior similar 
to experimental observation (Patterson [30]). Further, there are solutions (Shraiman 
k Bensimon [39]) that exist only for a finite time and culminate in a zero angled cusp 
at the interface. Howison [19] uses the class of known exact solutions to point out that 
the initial-value problem is ill-posed; it is possible to choose an initial condition for 
which the solution exists for all times, whereas there is a neighboring initial condition 
for which the interface develops a cusp after a finite time. 
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In essence, these theoretical results use a conformal map z((,t) which maps the 
interior of a unit semi-circle in the ( plane to the physical flow domain of a channel 
(see Figure 1). The location of the free surface at a time t is given by z(£, f) = 
x(9, t.) + iy{9 , t) for ( = e t0 on the arc of the semi-circle. The equation for the evolution 
of (x(9j),y(9j)) results from the usual application of the kinematic and dynamic 
conditions at the interface. The symmetry in the problem allows us to reflect the 
solution about the real axis. In essence, we include a mirror image of the channel 
so that the interface may be considered periodic. Thus, the conformal map must be 
analytic inside the unit disk, aside from a logarithmic singularity at ( = 0, but it 
may have singularities and zeros outside it. These can move towards the boundary 
of the unit disk, |£| = L In particular, zeros may reach the boundary in finite time, 
causing a cusp to form on the interface. The origin of ill-posedness, then, is that small 
perturbations can introduce a zero or a singularity in near the unit disk, which 
subsequently moves towards it and reaches it quickly. Following the work of Richardson 
[35] and Lacey [25], Tanveer [43] showed that all singularities, no matter what type, 
will move towards the unit disk, while preserving their type. Indeed, we surmise 
that round-off error in traditional numerical calculations introduces singularities at 
some distance from |£| = 1; these subsequently approach the unit disk and lead to 
the random pattern of tip-splitting, side-branching and finger competition seen in 
computations for increasingly small 5 . 

Ideally, one would like to understand the Hele-Shaw dynamics for a small non-zero 
B perturbatively, thus exploiting the simpler B — 0 case. There are several hurdles 
in accomplishing this. One is the ill-posedness of the 5 = 0 problem in the physical 
domain |(,'| < 1. Another is the fact that information in the 5 = 0 problem is not 
complete. We do not have an exact solution to the general initial value problem, nor 
do we have a detailed understanding of the motion of singularities. Furthermore, it 
is not known how to numerically compute solutions to the 5 = 0 problem effectively. 
Conventional numerical simulation in the physical flow domain \(\ < 1 suffers from 
uncontrolled growth of round-off errors (see Aitchison & Howison [1]). Employing a 
filtering procedure such as that used by Krasny [24] in the Kelvin-Helmholtz instability 
(another ill-posed problem) allows simulations in the physical domain to proceed. Still, 
as is demonstrated by Dai & Shelley [8], the choice of parameterization has a strong 
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effect on the accuracy of the computation. Moreover, a good parameterization is 
very dependent upon the initial data. For data with a general collection of zeros and 
singularities, reliable long time calculations in the physical domain become extremely 
difficult. 

In a recent development, Tanveer [43] has been able to demonstrate that analyti- 
cally extending the initial value problem for z((, t) into the region exterior to the unit 
disk leads to a well-posed evolution problem. It is important to note that when the 
initial interface location is known only to a finite precision, initial data 2((,0) in the 
extended region \{\ > 1 cannot be obtained in a well-posed manner. In effect, the 
ill-posedness of the dynamics is transfered to the ill-posedness of extending the data 
from |C| - 1 to |d > 1. Nevertheless, when data is specified in \(\ > 1 (say, with 2((,\0) 
given in closed form) the interface evolves without sensitivity to initial conditions. 

Tanveer's observation addresses the first hurdle mentioned above. In this paper, 
we address the second hurdle by presenting a numerical method which efficiently solves 
the initial value problem in the expanded domain when B = 0. Our work extends the 
method, developed by Baker & Tanveer [4], to include the trajectories of singularities 
of the form 

~ Mt) - Ut)) a * * 7*0,1, 2,... (1) 

explicitly in the complex plane ( when B = 0. Thus we are able to assess directly 
the impact of the close approach of pole and branch point singularities to |(| = 1. 
We do find that singularities induce tip-splitting and side-branching, and that the 
relative strength of the singularities control finger competition. Moreover, our studies 
draw connections between the parameters specifying the singularities and the result- 
ing physical behavior. In other words, for given initial data we are able to predict 
the outcome from knowledge about the initial singularities in \(\ > 1. Within this 
framework, comparisons with experimently observed features are possible by studying 
a random ensemble of initial conditions in \(\ > 1, subject to the constraint that they 
describe the same initial interface, up to some ^experimental' error. 

Beside the explicit treatment of singularities, there is another crucial ingredient in 
our method. An analytic function, such as a conformal map, is determined completely 
by its values on a closed curve. Instead of advancing the conformal map on |£| = 1 in 
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time as in standard boundary integral methods, our method advances the conformal 
map on a much larger circle. The interpolation of the conformal map to |(| = 1 is 
then a well-posed operation. In essence, knowledge of the conformal map on a much 
larger circle corresponds to knowledge of the conformal map on \(\ = 1 to a higher 
degree of precision. In contrast, extrapolation of the conformal map to a larger circle 
is an ill-posed operation, an operation not to be attempted numerically. Since our 
method uses the discrete Fourier transform to evaluate Laurent series, our results are 
spectrally accurate and can be performed in 0(N log N) operations. Furthermore, 
explicit treatment of the nearest singularities helps avoid point crowding typical of 
conformal maps. 

There are several limitations to our method. First, there is the matter of the initial 
data. We require that all singularities of in the extended domain be isolated branch 
points or poles. This rules out some types of initial data; for example, structures 
such as natural boundaries or essential singularities are not allowed in the extended 
domain. Nevertheless, the kinds of data we do allow produce a wide range of interfacial 
features such as tip-splitting, side-branching and finger competition which are similar 
to experimental observations. 

The second limitation is that, at present, our method is restricted to the zero 
surface tension problem. Given this, it is important to discuss the relationship between 
5 = 0 solutions and those for 5 > 0. The well-posed formulation of the 5 = 0 
problem allows the effects of small surface tension to be addressed perturbativelv. It 
turns out that, at the initial time, surface tension is a regular perturbation except in 
the neighborhood of zeros and certain singularities of Z(. Tanveer [43] has examined 
the effect of small surface tension on isolated singularities in the initial data of the 
form (1).* Capillary effects on singularities with a < — 4/3 cause only a regular 
perturbation. Their powers are unchanged but their speed and strengths A(t) are 
modified slightly by terms proportional to 5. Singularities with -4/3 < a < -1/2 
are immediately transformed into a localized cluster of -4/3 branch point singularities, 
though the behavior ( 1 ) is still relevant in an outer asymptotic sense. Thus, for times 
much less than 1/5 the interface behaves as though it is unaffected by surface tension, 

*Our use of o for the power corresponds to Tanveer’s [43] choice of — 0 . 
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provided these singularities do not come within a distance 0(B*+°) from |(,'| = 1 and 
zeros in z c are far from the unit disk. 8 For initial conditions with zeros present, recent 
evidence (Siegel, Tan veer, & Dai [40]) shows that small surface tension effects can 
significantly perturb the interface in 0(1) time. This can happen even when the zeros 
do not impinge on the unit disk. The surface tension effects occur in predictable ways 
when a localized cluster of -4/3 singularities created out of an initial zero (termed 
daughter singularity by Tanveer [43]) come within an 0(5 1/r3 ) neighborhood of the 
unit disk. However, there are initial conditions for which daughter singularity effects 
do not occur in 0(1) time. Our computations determine the leading order interfacial 
shapes in such cases. Perhaps more importantly, the method presented here provides 
a means to determine if and when the zero surface tension evolution of the interface 
deviates from the small surface tension evolution. 

In the next section, we describe the equations upon which our method is based. 
The explicit treatment of singularities is discussed in Section 3. Then we describe the 
numerical method in Section 4. In Section 5, we present tests of our method, and 
some typical results. Our conclusions are discussed in Section 6. 

2 The Equations of Motion 

In this section we present the equations which govern interfacial flow in a Hele-Shaw 
cell without the details of their derivation. We will follow closely the formulation in 
Tanveer [43]. Our discussion will be limited to flow in the channel geometry. For the 
equivalent formulation in a radial geometry, see Tanveer [43], 

Consider a Hele-Shaw cell of infinite length and finite width a in which air of 
negligible viscosity is pushing a viscous, incompressible liquid. Introduce the conformal 
map a(C, t) which takes the interior of a unit semi-circle in the ( plane into the viscous 
fluid region of the channel, which lies in the z plane. The circular arc |C| = 1 is mapped 
to the interface, and the diameter is mapped to the channel walls. A schematic of the 
mapping is provided in Figure 1. Note that we set the width a = 2. 

5 When singularities are initially Of 1 ) distance from the unit disk, they can move to within 0( B a+S' ) 
in 0(ln ) time. 
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The functional form of the conformal map is given by 


2 

*(C*0= In < + * + /(C,J) (2) 

7T 

where / is analytic in an open set which contains the unit semi-circle. The analyticity 
of / on the circular arc guarantees the smoothness of the interface. We require that 


lm{/} = 0 (3) 

on the real diameter of the semi-circle to satisfy the condition that 2 maps this diameter 
to the channel walls. In addition, we assume that ^ 0 in a region containing the 
unit semi-circle. The Schwartz reflection principle then implies that / is analytic and 
2^/0 for |C| < 1. 

The fluid velocity u, averaged across the plate gap, satisfies Darcy’s law 


u = 


-b 2 

12 p 


Vp 


where p is the viscosity, b is the plate gap, and p is the pressure (here considered 
as a function of x and y). Thus (-b 2 /V2p)p provides a velocity potential <f>. Incom- 
pressibility implies the existence of a stream function ip. We can therefore introduce 
a complex potential function W(z,t) = 4>(z, t) + iip(z,t) which is an analytic function 
of 2 in the fluid region of the channel. Considered as a function of (, this potential is 
decomposed as 

W ( Ci t ) = — ~ In £ + i + u( (, , t ) ( 4 ) 

where u> is assumed to be analytic in unit semi-circle, implying (via the Schwartz 
reflection principle) its analyticity for |(| < 1. In (4) the velocity at infinity is assumed 
to be 1; together with a = 2, this choice makes our variables effectively dimensionless. 
The relation 


Im{u;} = 0 


(5) 


is required to hold on the real diameter of ,9, and is the mathematical statement of 
the condition that there is no flow through the walls. 

The interfacial conditions will determine the evolution of the map z((.t). The kine- 
matic condition states that the normal component of the fluid velocity is continuous 



across the interface and implies that on \(\ — 1 


Re 




( 6 ) 


Details of the derivation of this equation are available in Saffman [37] and Richardson 
[35]. In the absence of capillary effects, the dynamic condition requires continuity of 
pressure across the interface. When the viscosity of the less viscous fluid is negligible, 
this condition gives 

Re {a;} = 0 (7) 


on |(| = 1. We note in passing that a more complete description of Hele-Shaw flow 
introduces more complicated interfacial conditions due to a thin film left behind by an 
advancing interface. However, several numerical (Park & Homsy [29], Rein el t [33, 34]) 
and analytic (Tanveer [44]) studies have shown that many salient features of interfacial 
motion described by the more detailed model are captured by the simpler interfacial 
conditions (6-7). 

From (5) and (7) it is immediately apparent that u; = 0 for all ( in the complex 
plane. Thus, (6) becomes 


Re 


( 8 ) 


1 = 2_ 

'kl 2 

In order to extend (8) into the region |(j| > 1, we first provide an analytic extension 
of the conjugate of an analytic function evaluated on | ( |= 1. Let 


no = J2 

n~ 0 

be an analytic function in |£| < 1. Then 

oo 

/(C) = £«*(,*” (9) 

71=0 

is analytic in |(| < 1, where the overbar denotes complex conjugation. For functions 
that are not necessarily analytic at the origin, though analytic on a segment of |(,j = 1, 
we can generalize the above definition of F(Q through the relation F(() = /’(C). 
/"(I/O is then the analytic extension of F from | (j |= 1. It follows that 


Ke{F(C)} = ^[F(C) + /(l/C)] 


( 10 ) 



for ( = e' e . In addition, |F(()| 2 = F(QF(l/0 on (,' = e ’ e . Thus (8) can be written as 


Re 


Jt_ 

( z < 


(l/d*) 


ID 


on ( = e' e . 


The continuation of equation (11) into the domain |d < 1 can now be obtained in 
a straight forward manner bv employing the Poisson Integral Formula. In particular, 
we use a variant of the standard formula which gives the value of an analytic function 
in the domain |d < 1 in terms of its real part evaluated on the unit circle. Application 
of the formula to the function Zt/(C z O Y^Ws (with appropriate choice of imaginary 
constant) 

( 12 ) 


7^ = /«,V) 


for |(| < 1, where 


7T ‘ J\C\=l 


1 


C + CdC 


13) 


*<(CM)*<(iAV)C'-C C' 

Equation (12) can be analytically continued into the domain |(| > 1 by deforming 
the contour in the usual way, producing an additional term from the residue of the 
pole at (. Consequently, we have 

2C 


z t = C - 


(14) 


i/C»<) 

for |d > 1. Note that /((>!) defines different analytic functions in |(| < 1 and |(| > 1. 
A useful alternative form for ( 14) may be obtained by employing (10) to write (8) 


as 


C z < 


Qztd/CJ) 2C 


15) 


Since z t (l/Q and ^(l/() are analytic in | ( |> 1, this equation provides us with a 
second expression for /((,/), 


/«•„=_< Mim 

when |(| > 1. 

In order to make the structure of (12) and (14) more apparent, define 


16) 


91 

92 


C/(dO 

— 2( 

5 c (i/0* 


(17) 

(18) 



Then (12) becomes 


(19) 


*t ~ q\zc = 0 . 

for |C | < 1, and (14) becomes 

zt~q i*< = ft- (20) 

for |C| > 1. These equations, corresponding to equations (3.1) and (3.3) of Tanveer 
[43], have the advantage over (8) of allowing studies of the presence and influence of 
singularities of in |C| > 1. Furthermore, they lead to the development of a new 
numerical procedure for calculating interfacial motion in a well-posed fashion for a 
restricted class of initial conditions. 

Since /((,*) and z c (l/C) are analytic in |C| > 1, so too are <7i(C,<) and q 2 ((,t), 
except at infinity, where q\ has a simple pole, i.e.. grows as Q. Since (20) has a 
form analogous to a first order hyperbolic system in the complex plane, (although in 
reality it is a nonlinear integro-differential equation with coefficients q\ and q 2 that 
depend nonlocally on r), the analyticity of qi and q 2 has important consequences on 
the presence and motion of singularities of 2 . For the convenience of the reader, we 
provide a summary of what is known (Richardson [35], Lacey [25], Tanveer [43]): 


1. There is no spontaneous generation of singularities in the finite complex plane. 
Furthermore, the form of a singularity which is present initially in the region 
Id > 1 is invariant with time. Singularities which are present initially at infinity 
do not move to a finite location. 


2. VVe define a “characteristic” in the ( plane by 


rfCc(0 

dt 




( 21 ) 


Let d( t ) = RA t) then 




-Re 


f ¥i(Cc(0,0 ] 
\ C c«) J 


( 22 ) 


Tanveer [43] has shown that the right hand side of (22) is less than zero when 
Id | > 1. Consequently, “information” outside of the unit disc flows inward 
toward |d = 1. 
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3. This is particularly true for isolated singularities of the form (1). The location of 
these singularities satisfies (21), i.e. they move with speed -<7i(Cs(*Mh so they 
move towards the unit circle. Incidentally, the singularity in z which is present 
initially at ( = 0 does not move in time, since the characteristic speed at ( = 0 
(given by <71 (0,0) * 1S zero. 

4. Singularities with a > —1/2 reach the unit circle in finite time. Singularities 
with a < —1/2 come indefinitely close to, but never reach the boundary |C| = 1. 

Several properties of zeros in are also relevant. For example, there is no spon- 
taneous generation of zeros of z<^ in the finite complex plane, although zeros which 
are present initially at infinity can move to a finite ( location. When a zero impinges 
on |d = 1, it produces a zero-angled cusp in the shape of the interface; see Howison 
[19. 20] and Shraiman & Bensimon [39] for exact solutions where this happens. There 
is no physically sensible way of continuing the solution in time. Unlike singularities, 
the zeros of z< do not move generally with the characteristic speed -q\ and it is hard 
to predict if a given zero will hit the unit circle or stall at a finite distance from |d = 1. 

In this paper, our interest will be focussed on initial conditions containing only 
isolated singularities of form (1) with a < -1/2, so they do not reach |d = 1 in finite 
time. We shall also pick cases where zeros in z^ do not reach |d = 1 during the time 
of our computations. For these cases, the inclusion of surface tension acts in a regular 
perturbative manner, so the results indicate what can be expected in the limit of weak 
surface tension. 

3 Explicit Treatment of Singularities 

Baker k Tanveer [4] use (20) to solve directly for z((, t) on a circle in the ( plane of 
radius R(t), assuming that this circle does not contain a zero of z^ or a singularity of 
- z-i + (2/tt) In C- The advantage of computing the evolution of f((J) on the 
boundary |(| = R(t) is that during the process no singularities in z<; will be introduced 
in |Cl <#(£). Of course, singularities in z^ may be present outside |£| = ft(f), but 
as long as R(t) shrinks fast enough, singularities will not intrude into |£| < R(t): the 
presence of zeros must be checked separately. Furt hermore, a function which is analytic 
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in |C| < R(t) can be determined throughout this region from its values on |£| = #(/) 
through evaluation of its Taylor series - this step is well-posed as demonstrated in the 
next section. Thus, the interface, z(|C| = Id), can be recovered from the solution on 

Id = R(t)- 

Unfortunately, R(t) must shrink faster than the approach of any singularity to 
Id = 1. Consequently, the computation may terminate (when R = 1) before any 
interesting structures have developed on the interface. For initial conditions containing 
isolated singularities of form ( 1 ), this obstacle may be removed by explicit treatment of 
the singularities. In this section, we derive equations in which singularities in |(| > 1 
are treated explicitly; they are in a sense ‘subtracted from the data and evolved 
separately. In the next section, we show how the resulting equations can be solved in 
an efficient and well-posed fashion. 

For initial data containing singularities of form (1) in |d > C the solution may be 
written in the form. 


J i 


Z (CJ) = 

j= i ' 

+ #«,<)( i 


. \ «,+i j 

1 ' + E 

i +i 


c 


WJ 

\ °J +1 


a, + 1 


U-oio, 




+ G((,t) InC 

7 r 


(23) 


where |C ; I > 1; (*j are real constants (excluding 0,1,2,...); and J ( Cj ( 0 ) , 0 ) / 0. 
Here (j(t) for j = 1, . . J\ mark the location of singularities on the real line, and (j(t) 

for j = J\ + 1 J indicate singularities with non-zero imaginary part. In order to 

satisfy (3), the amplitudes E j (£,t) for j = 1 h must be real. The amplitudes 

EHX-t) for j = J\ + 1 J are, in general, complex. Conjugate singularities Cj(f) 

with amplitudes are included to satisfy condition (3). The branch cut is 

chosen so that the argument of (1 — (,/(,j) li p s between ?r and —i r. we have chosen 
this form to simplify the numerical evaluation of the interfacial location when | C I = 1- 
In the region |(| > 1, there can be other singularities at (j which are not explicitly- 
represented. However, it is assumed that the are the nearest singularities so that, 
given i, the inequality |(,',j > R(t) > | Cj I holds for all j = 1,. The functions 

E J ((,t) and G'U'.O are analytic in the annulus 1 < |C| < R(1)- If all singularities are 
explicitly represented, these functions are analytic in the entire region exterior to the 
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unit disk. Unfortunately, it is possible for E J and G to have singularities in \(\ < i 
which cancel in (23) so that z + ( 2 / 7r ) 1 n C is analytic in \(\ < 1. 

In order for (23) to be a solution to (20), (j, E J , and G must satisfy 


dQ 


dt 




(24) 


E J t ~ 9i 


E\ = (oj + 1) E 3 { (25) 

c 3 \ C-Cj(t) 0(0 J 


Gt - <7i G( = — — -q\ + q 2 - 


(26) 


In the case of initial data with logarithmic singularities (poles in Z ( ) the solution 
has a modified form, 


*(C»0 = j>(C,0 In 1- 


c 


j = 1 

+ E J «;j ) In ( 1 


0 ( 0 , 


0 ( 0 , 


J 

+ E 

j=Ji+l L 
2 


E 3 ( 00 In 1 


<, 


0 ( 0 , 


where 


^0 

£0 


+ 0'(4, 0 In (, 

x 


-<?i(0(o,o 

9i E\ = 0 


Gt - 9l G( - qi + ]T] 

j=! 


<71(00 - «7i(00)O) 9l(0O) 


(27) 

(28) 

(29) 

(30) 


C - 0(0 0(0 J 

Of course, it is possible to have a combination of the two forms, (23) and (27). 

Our numerical method is based on (24-26) or (28-30): in the next section, we 
describe how the solution may be updated in a well-posed fashion. 


4 Numerical Method 

Our numerical treatment is a combination of tracking the positions of the singular- 
ities by solving (24) (or (28)) and of advancing E J ((,t) and G((.t) in time by the 
method of lines. Despite appearances, it is not necessary to update and 

G((,t) throughout the computational domain, an annulus 1 < \(\ < R(t). Knowl- 
edge of these functions on the circles \(\ — 1 and |£| = R(t) is sufficient to determine 
them everywhere inside the annulus. Nevertheless, we must use a special procedure 
to update E J ({, t) and G(C>0 in a numerically stable way. 
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We introduce a decomposition for a function #(<) which is analytic in the annulus 
by writing it as the sum of a function g+(0 which is analytic inside a circle |C| = ft 
and a function g~{0 which is analytic outside the circle |C| = 1- The decomposition is 
made unique by requiring ff+(0) = 0. By considering the Laurent expansion for g( Q, 
we note that 

g-(0 = H ( 3 n 

fc= — oc 

g+(0 = ^2$kC k < 32 > 

i 

where the Taylor series for £+(<,*) will have a radius of convergence ft + > R , and the 
Laurent series for g~(0 will converge outside |C| = R- < T We call the function 
g + inner analytic and the function g„ outer analytic. 

An important feat ure of an inner analytic function is that its values inside a closed 
curve can be determined from the values on the curve in a well-posed fashion. To see 
this, consider an inner analytic function g+(Q evaluated on a circle of radius R. Let 

(,* = Re ' 8 , then 

g + (9) = f^g k R k e' ke 

k= 1 

In practice, the coefficients g k R k are obtained by the discrete Fourier Transform. Since 
the sum converges, these coefficients must decay with increasing k. There will be a 
value, k r say, where the coefficients reach the round-off levels of a computer. If suffi- 
cient resolution is used with the discrete Fourier Transform, then the coefficients with 
k > k r will contain round-off levels and not their actual values. This has important 
consequences when we wish to evaluate g+ on a different circle Q — re . Then 

so each coefficient g k R k is multiplied by r k / R k prior to using the discrete Fourier 
Transform again to evaluate the sum. If r > ft, then round-off errors in g k R k for k > k T 
will be significantly amplified, even to the point that all accuracy in g + (C = re' 9 ) will 
be lost. However, if r < ft, round-off errors are decreased. For an outer analytic 
function the situation is reversed: evaluation on r > ft (here we will take ft = 1 ) is 
numerically stable, while on r < ft it is not. The exception to our statements is when 
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there is a known, finite number of terms in (31) or (32) with amplitudes well above 
round-off levels. Then the relative errors in the evaluation of g+ or on any circle 
remains small when we use only the finite number of terms. 


4.1 Basic Algorithm 

We decompose both £V((,\<) and G((,t) into inner and outer analytic parts and use 
the method of lines to update E J + , G+ on |(| = R , and E J _, G_ on |C| = 1. We 
obtain evolution equations for these quantities by applying the projection operators, 
H+f = /+ and H-f = /_. Thus, from (25) and (26), we have 


E 


: 


H. {/?,«,<)} 

H- i 91^ + (®j + 1) 


f 9i(cV) - 9 i(Q(<M) 

V C-0 


9i (Cj,p 
0(0 



(33) 

(34) 


6'_, — — ^9i(0 0 + 92(<o 0 j • (35) 

Before giving the results for the other components, we note that q\ and <72 have very 
simple decompositions. From (16) and (17): 


91 (00 = C<7i(0 + 91- 


(36) 


where q\ + = (q\(t) contains only one term, whereas ( 18) shows that <72 is outer analytic 
only, 

92(0 0 — ?2_ • (37) 

The form of the decompositions (36) and (37) is very beneficial in our design of a 
well-posed algorithm. For example, we use the facts that the product of E 3 _ with (71 
and the product of E J _ with [91(C) - 9i(0t)]/(C - Ofc ) result in only outer analytic parts 
to simplify the H+ projection of (25): 


E ] +t = W + {i? 2 (00} 


<, ~ 0 


0(0 


The resulting equation for G+ simplifies to 

G+ , = 'H+ {91^+,;} 

when we use the facts that -2<7 i/(7tC), 92 and q\G-. are outer analytic. 


(38) 


(39) 
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The following properties can be easily deduced from (34-35) and (38-39): 

1. The evolution of E+ does not depend explicitly on E J _. As a result, if £*>(£, 0) 
satisfies £^.(£, 0) = 0, then E 3 + (£,t) = 0 for all t. The same result holds for 

G+(C,<). 

2. If £^(C,0) contains a finite number of modes in (32) with highest wavenumber 
k — M , then no modes with wavenumber k > M will be generated in E+(£,t.). 
The same result holds for G+. 

Similar results hold for the decomposition of E 3 and G in (29) and (30) for loga- 
rithmic singularities. 

To apply the method of lines to (34-35) and (38-39), we assume that we know E+ 
and G+ at N evenly spaced points on the circle r = R(t), and E J _ and 6 f _ at TV evenly 
spaced points on the unit circle. Actually, it is not necessary to use the same number 
of points on both circles, but we make that assumption for ease of presentation. First, 
we describe how to evaluate the right hand side of (34) and (35) on the unit circle. For 
the moment we assume that q\ and q 2 are known there; we describe their computation 
in detail in the next sub-section. 

We compute the coefficients CkR k in the representation 

N/2 

E } + (<; = Re' e ,t) = Y t c k R k e lk9 

k=l 

by use of the discrete Fourier Transform. The coefficients c are obtained by a division 
by R k and then used to evaluate £+ and E+ on the unit circle by use of the discrete 
Fourier Transform, where 

N/2 

£+<(c,o =£> c *c* -1 - (40) 

k= 1 

The function E J _, is also determined on |(| = 1 by using the discrete fourier Transform. 
Then we have all the information needed to evaluate R\ in (34) at N evenly spaced 
points on the unit circle. To execute the projection W_, we use the discrete Fourier 
Transform to calculate the coefficients in the representation 

N/2 

Ri(C<)= E r*UK fc . 

k=-N/2 


16 



Then we set to zero all the coefficients with k > 1. The result is an outer analytic 
function, which is the forcing term in (34). Balance of the Fourier coefficients of like 
modes on the left and right hand sides of (34) then yields a set of equations 

-^(0 = r k (t) 


for the Fourier coefficients E J k (t) ( k < 0) of E 3 _. In an equivalent procedure, we may 


use G'+, known on r = R{t), to determine G + and use G'_, known on the unit circle, 
to get G-^. Then the forcing term for (35) is computed by executing the projection 
H- as described above. 

From our discussion on the importance of “characteristics” defined by (21), we 
anticipate that E J + and G+ must be evaluated on a circle with a radius that is collapsing 
at rate. 


1 dR 

-E~T7 = - max Re 
R dt |<|=R(<) 


f gi(C,0 j 


(41 


Then, information outside the circle |(j| = R(t) will not cross the boundary into the 
annulus, 1 < |(j < R(t). We include an advection term that accounts for the change 
in E J + due to the change in R(t): thus (38) becomes 


dE 3 + 

dt 


= H + 


#2(C. t) + 


^dR 

H~dt 


T} 


( 42 ) 


A similar term is needed in the equation (39) for G + . 

We will describe in detail in the next sub-section how qq and </ 2 may be evaluated 
on (C| = R{t). The quantities E+ , G +< can be evaluated by spectral techniques as 
described in (40). The expression contained within brackets on the right hand side 
of (42) can then be evaluated. The Fourier coefficients of the right hand side of (42) 
are obtained using the discrete Fourier Transform, and the projection H+ is done by 
zeroing all modes with k < 0. Equating like Fourier modes then leads to a set of 
evolution equations for E J k (t), the Fourier coefficients of E } + , for k > 0. The same 
procedure may be used on the equation for G+. 

Any suitable ODE Solver may be applied to the evolution equations for Q, El, and 
Gk ■ For the results used in this paper, we use the standard fourth-order Runge-Kutta 
method with fixed step-size. 

Since 2(<j, 0 + (2 /tt )ln C has a Taylor series expansion about ( = 0, all the negative 
terms in the Laurent expansion for G((j,<) must cancel all the negative terms in the 
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Laurent expansions of the sums in (23). This will be only approximately true in the 
discrete calculation, so it can be used as a check for numerical accuracy. Although it 
is not necessary to compute to update the interface, we have done so in order 

to employ this check. 

For certain initial conditions, our algorithm is greatly simplified. For example, 
when E+ has only a few terms in its Taylor series expansion, then it is possible to 
write down by hand the evolution equations for the Fourier coefficients, E 3 k , since the 
products 


^ f (lR \ ri 
qi + Rli) E+ < 


and 


<7iU\0 - 9i(C. v f ) 




C -0(0 

will contribute only a few terms due to the simple decomposition (36) for q \ . In 
particular, if £+, E i, G+, and are time-dependent constants in then only E 3 _ 
need be computed. 


4.2 Computation of q\ and q 2 


Here we provide details of the computation of q\ and q 2 in \(\ > 1 bv using (13), (17), 
and (18). First, we calculate £|((, f ), and E^^CJ) on \(\ = 1 by means of 

the discrete Fourier Transform as described in the previous subsection. These values 
are used to compute 1/^(0 — PooiO/ D(() on \Q — 1 where 


D( 0 = - 


2ppo(C) 


J 1 


- Y^pmo 




(«> + !) 


m o 


- E < (0 1 - i 


poM 'i 

j=Jl +1 

+ C<(0 Poo(C) 






wj 


(43) 


and 

\ -a i J 

^)=nh n 

,=i v c/ ,*j 1+1 

1^3 

This form for l/z((C) is obtained by differentiating (23) with respect to (,\ then fac- 
toring out l/poo(0 and taking the reciprocal. We have found it necessary to compute 
1/^((J) using this expression to obtain accurate results for q\ when singularities are 
very close to |£| = 1. The function l/z((\/Q can be computed on \£\ = 1 by taking 
the conjugate of 1/^(0- Consequently, the function q 2 can be computed on the unit 
circle using ( 18). 
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The function q\ is computed as follows. We write 


2 


= do + ^2(dk( k + dk rC k ) 

k=i 


( 44 ) 


on ( = e' e . We determine approximations to the N coefficients 
jd] , . . . , d/v/ 2 ) by means of the discrete Fourier Transform. Upon 
(44) into ( 13) and computation of the residues, one finds that 

f N/2 1 


9l = ~C 


do + 2 ^ d k C k 


k= 1 


. . . , d A y 2 }^ 

substitution of 


(45) 


Note that d 3 is real (i.e., d 3 ~ dj ) for the channel geometry, as a consequence of z being 
real along the real axis of (. The functions q\ and <72 can be analytically continued out 
to the circle |C| = R(t) by the ideas expressed in (40). Similarly, q x can be evaluated 
at the location (j(t) of the singularities. 


4.3 Computations for Extremely Close Singularities 

Many of the interesting interfacial features revealed by our calculations occur when 
singularities approach the unit circle within a distance much less than machine preci- 
sion. To track these singularities reliably when they are that close to the unit circle 
we must express their location as 

0 - ( 1 + y ej 

The time evolution of the quantities 0j(t) and A 3 (t) - 1 jh 3 are then determined 

We use A j instead of 6j so that singularities can be allowed to come extremely close 
without the worry of them jumping inside the unit circle due to time-stepping errors. 

The computation of q\ ( (jj , /) is more delicate when (j is extremely close to the unit 
circle, since the quantities 

1 _ 1 

FcTruT^ 


numerically using the equations 

dA, 

dt 

dt 
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in (45) will be approximated as 1, leading to inaccuracies in the computation. In order 
to obtain a better value for 1/Cf we write 


1 

1 + b j 


= 1 


h 

1 + b j 


1 + /l ( ) 


and compute fk(bj). defined by 


1 

(1 + «i)* 


1 + fk{bj ) , 


(46) 


using the recursion 


/*(*;) = 


1 -f- 6j 


The substitution of (46) into (45) then yields = qi(e' e > ,t) + 6qi(Cj,t), where 


N/2 


bq^CjJ) = -e ie > I bjd 0 + 2^4 [6 } (1 +/*) + Ik] e~ ,k0J 


(47 


k = 1 


Then, 





(48) 


where we have used 

I ; J ^(O^O/C) 

when |C| = 1, which follows from (13) and (17). By using (48), we avoid any problems 
with round-off errors. On the other hand, Im{gi/d} is computed with q x directly 
determined from (47). 


5 Numerical Results 


We use a known exact solution due to Saffman [37] to help validate our method: 

2 i ( e x 

z(Q.t) = i + d(t) - - In ( + ~ In I 1 - 

7T 7T V 


a 2 (t), 


where 1 < a(0) < oo. The functions d(t) and a(t) are determined by 

, 1/4 


a(t) = \ + K\e 


-2rr t 


(49) 


and 


d( t) = A o + 2f + 7T“ In [l + A i e 

In L 


-2 TTt 
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where A'o and A*i are constants which are determined by the initial conditions. From 
(49), it is clear that has poles at Ci ,2 — ±a(/), but no zeros. The form of the initial 
data (27) used in our numerical method is made to correspond to (49) by setting 
E 1 ( C,0) = E 2 ( C,0) = 1/tt and G'(C, 0) = i + d( 0). We pick d(0) so that the initial 
profile has zero mean height. The functions E J ((,t) and G{(,t) are constants in ( for 
this problem, and therefore only a single Fourier mode (namely the constant mode) 
is necessary to specify them. Nevertheless, we use 512 modes in our computation in 
order to check the numerical stability of our algorithm for computing E J (£,t) and 

G((, t). 

The results of our numerical method for a(0) = 2.0 are shown in Figure 2. We use 
a time step of At = 0.005. At / = 3, the calculation gives the positions of the poles as 
Ci , 2 = ±1.0000000244216, whereas the exact positions are Ci,2 = ±1.0000000244215. 
Although the difference is less than 10~ 12 , the exact and the numerically computed 
profiles differ by only a little less than 10 -6 . However, we find no growth of round-off 
errors in the modes comprising E J ((,t) and In this example we set R( 0) = 

2,000, and stop the calculation at t = 3 when R is nearly 1. Because there are no 
zeros in we can take A(0) much larger and run for even longer times. 

There are no known exact solutions for channel geometry when branch point sin- 
gularities are located in the region |(| > 1 (although exact solutions with branch point 
singularities have been found for sector geometry; see Tu [47]). We replace the two 
pole singularities in of (49) by two branch point singularities of equal power and 
amplitude, and locate them symmetrically on the real axis of C- Our choice is mo- 
tivated by the desire to understand the role of the power of a singularity in z ^ on 
the shape of the interface. Thus we select initial data corresponding to (23) with 
J\ = 2, Q| = » 2 , and (\(0) = — C2(0). In Figure 3, a\ = —4/5, and in Figure 4, 
ai = -4/3. Initial values for the singularity positions £j(0) and amplitudes £ J '(£,0) 
are determined through experimentation; those that lead to solutions in which the 
interface becomes well-deformed before zeros impinge on the unit disc are selected 
for presentation. In this manner, the initial singularity positions and amplitudes 
are chosen as Ci ( 0 ) = — C 2 < 0) = 1.2, E l ((,0) - £ 2 (£,0) - 1.94 for Figure 3, and 
Ci(0) = — C 2 ( 0 ) = 3.3, A J ((,0) = £ 2 (C,0) - -2.2 for Figure 4. In these and all 
subsequent computations w r e set G((,\0) = * + c, where c is a real constant selected to 
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produce an initial profile with zero mean height. Recall that for constant initial data 
such as this, only negative Fourier modes in E } and G are generated. Consequently, 
we only need to solve equations (34) and (35). The positive Fourier modes are set to 
zero after each time step in the calculation. We use 512 points on the unit circle in 
the plane, and a time step of A t = 0.0005. 

We check our numerical results by comparing them with the results obtained when 
none of the singularities outside the unit disk are represented explicitly. To perform 
the latter calculation, we use the code to solve for G'+(C, t) + Go(t) on a circle of radius 
R(t) that is closer than the smallest |(,j(/)|. Consequently, we shrink the radius R(t) 
as described in section 4. While the calculation must be stopped before long because 
R < 1, the output provides a useful check on the present results for some period of 
time. In addition, we decrease the time step and increase the number of modes until 
there are no detectable differences in the solution within plotting accuracy. 

Unlike the previous example, we cannot be sure that zeros do not approach the 
unit disk. So we monitor the presence of zeros by computing the integral 


N x 


-Li 

27h y|<|=i 



(50) 


which, according to the argument principle, equals the number of zeroes minus the 
number of poles of z ^ inside the unit disk. As long as there are no zeros of z<; in 
|C| < 1, this integral will equal -1, due to the simple pole of Z(_ at C = 0. For 
the results presented in Figures 3 and 4, the value of N z remains within 1% of -1 
throughout the length of the computation. Note further that if a zero in z <; is located 
near the unit disk, the close presence of a pole singularity in the integrand of (50) 
causes a loss of accuracy in its numerical evaluation. Thus we feel confident that no 
zeros are very near to the unit disk up until the times of the final profiles in Figures 


3,4, and 5. 


The contrast between Figures 3 and 4 is quite striking. The higher value of 
a = —4/5 in Figure 3 correspond to a singularity of weaker effect in that the bulge of 
fluid at the base of the finger is less rounded, causing a finger that is more pointed. 
The lower value of a - -4/3 in Figure 4 produces a more spherically shaped bulge of 
fluid, causing a thinner neck at the base of the finger. In general, our experiments with 
various powers a show that for powers a > -1, more pointed fingers are produced 
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from the less pronounced bulges of fluid at their bases. For a < -1, the fingers have 
parallel sides, but their bases have thinner necks as a consequence of more spherically 
rounded bulges of fluid there. 

A representative Fourier spectrum from the calculation with a = -4/3 is pre- 
sented in Figure 5. There is no sign of spurious growth in the large | Ar| modes or any 
other indications of numerical instability, although it is necessary to evaluate q\ using 
expression (43) in order to avoid such growth. The rise in the tail of the spectrum 
with respect to time is due to singularities in the functions E J , G , q\, and q 2 at po- 
sitions £ = 1 /Ci and £ = l/( 2 inside the unit disc. Since z - (2/jr)ln£ is analytic in 
|C| < 1, the singularities in E J and G must cancel out in the expression for z(£), (23). 
Nevertheless, these singularities do affect the numerical computation of the quantities 
E\ G, <71, or q 2 • Since the singularities move toward £ = ±1 from inside the unit 
disc as £1 1 2 — * ► ±1 from the outside, a large number of Fourier modes are required 
to obtain accurate representations for these quantities when the singularities are very 
close to the unit circle. Another consequence of the close approach of the singularities 
to the unit circle is the need for much smaller time steps to maintain accuracy. This 
is also true when zeros in z ^ get close to the unit disk. In general, it is either the 
close approach of singularities in E J and G or zeros in z^ to | £ | = 1 that force us to 
terminate our calculations. 

If we ignore the loss of accuracy, we can continue the calculations shown in Figures 
3 and 4 a little further in time. We see evidence that zeros are impinging on the unit 
disk. For the calculations associated with Figure 3, a zero is approaching |£| = 1 at 
a point corresponding to the tip of the finger. Thus, we expect the finger to form a 
cusp in finite time. For Figure 4, two zeros approach |£| = 1 at points corresponding 
to the tops of either side of the finger, so we expect cusps to form at these positions. 
Asymptotic theory suggests that the initial zeros will give rise to localized clusters 
of -4/3 singularities (daughter singularities) when 0 < B << 1. The leading order 
motion of each of these clusters satisfies equation (21), i.e., a cluster located at Q 
moves with speed -q\{Q{t)J). If such a cluster comes close to |£| = 1, it can cause 
the interface to deviate significantly from the B = 0 solution. 

We briefly consider the influence of non-zero surface tension on interfacial shapes 
by comparing the zero surface tension solution to that for small surface tension. The 
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non-zero 5 solution is obtained using the boundary integral method developed by Hou, 
Lowengrub & Shelley [17]. In the first case (Figure 6a), we consider initial data with 
zeros of initially placed at infinity, and with singularities which satisfy aj < —4/3.^ 
The value of surface tension is set to B = .0025. From asymptotic theory, it is expected 
that the addition of a small amount of surface tension will make little difference in the 
evolution of the interface for at least 0(ln B) time. The actual agreement between the 
interfacial shapes is quite remarkable: the two solutions are indistinguishable over the 
entire length of the run. The agreement is unaffected by using real singularities which 
only satisfy aj < -1/2, since these singularities behave as though they are essentially 
unaffected by capillary effects for the time of the computation. 

In contrast, when the initial data of Figure 4 is used, the difference between the 
5 = 0 and 5 = .0025 shapes is very significant. As seen in Figure 6b, the B > 0 
finger eventually diverges from the corresponding zero surface tension solution and 
approaches a broad, steadily propagating finger. The broadening is apparently caused 
by the approach of daughter singularities created from an initial zero in \(\ > 1 (see 
Siegel, Tanveer & Dai [40]). 

The above examples show that, in some cases, surface tension causes a singular 
perturbation in 0(1) time, whereas in other cases it does not. Our method gives 
us a means of discerning these cases through computation of daughter singularity 
trajectories in the extended domain [43]. We remark that our attempts to calculate 
the 5 = 0 solutions using the boundary integral method with filtering failed before 
the interface had advanced very far. This is because the close presence of strong 
singularities causes fast growth of the high wavenumber modes, and numerical noise 
quickly contaminates the computation. Thus, boundary integral methods often appear 
unsuitable for comparing 5 = 0 solutions to those for 5 > 0 over times in which the 
interface becomes significantly deformed. Dai & Shelley [8] report related problems in 
5 = 0 calculations, as discussed in the introduction. 

We turn now to a consideration of the influence of additional -4/3 singularities 
with weak amplitudes during finger formation. Figure 7a illustrates the interface 
evolution resulting from initial data of the form (23) with J\ =2, J = 3 and with 

^This is accomplished by prescribing initial data in z<; of the form 2<(C<0) = —2(1 - C 2 /C? P/f^O* 
Here we use a = — 3.2 and Cs = (1.6. 0.0). 
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singularity strengths Qj = -4/3 for j — 1,...,3. The initial singularity positions 
were chosen as Ci(t) = — C2G) = (1.65,0.0), and (3 {t) = (—0.6,1.83), and the initial 
amplitudes as E l (C 0) = £ 2 (C 0) = (-0.9, 0.0), and £ 3 (C,0) = (-0.01,0.03). The 
calculation uses ^ = 512 points and a time step At = 0.0005. Two of the singularities 
((1 and (2) reside on the real line, and these have large enough amplitude to produce 
the wide bulges of fluid centered at y = 1 and y — — 1 which define the main finger. 
As the third singularity approaches the unit disk, it generates a small bulge of fluid 
that gives the appearance of the formation of a dimple on the evolving finger. 

The dimple is clearly stationary in the laboratory frame, despite the overall growth 
of the main finger. Such behavior has been well documented in laboratory experiments 
( e.g. Park & Homsy [28], Tabeling, Zocchi, & Libchaber [41]) and in numerical 
calculations (e.g. DeGregoria and Schwartz [9, 10]) with B ^ 0. Our results show 
that it is the specific nature of the trajectory of the singularity £3 (t) that accounts 
for this behavior. We show the trajectories of all the singularities in Figure 7b. At 
first, a dimple starts to form as &(t) approaches the unit disk. As (3 ( f ) begins to 
move around the boundary of the unit disk towards ( = -1, the dimple continues to 
grows but remains stationary in the laboratory frame. With the assumption that the 
singularity is close to either £ = ±1, we can show that its speed is actually the correct 
speed for it to remain stationary in the laboratory frame. Details will be provided 
elsewhere. The point to be made here is only that certain physical properties can 
be understood in terms of the motion of singularities in the complex plane. Under 
the presence of surface tension, the narrow finger will eventually widen; however, the 
formation of the dimple and its relationship to the motion of the singularities will not 
be affected. 

Placing additional singularities in the complex plane produces additional dimples 
and leads to the appearance of side-branching. An example of side-branching due to 
multiple pole singularities is given in Figure 8a. We use pole singularities, rather than 
a = —4/3 singularities, since our ability to track them arbitrarily close to the unit 
disc leads to a more dramatic example of sidebranching, but in reality side-branching 
is more likely to occur from the patterns of —4/3 singularities created by the initial 
transformation of zeros in z^. Initial data corresponding to (27) with J\ = 2 and J ~ 14 
is used to generate the profiles. The starting amplitudes E J ((, 0) and singularity 
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positions Cj(0) are given in Table I. The two main singularities and ( 2 produce the 
finger centered in the channel, and remaining singularities Q for j = 3, J of smaller 
amplitude are located to cause side-branches to form near the base of the finger. 
Additional singularities can be placed to allow side-branching to run the length of the 
finger. As expected, the dimples and the corresponding side-branches are stationary 
in the laboratory frame. The motion of the singularities is shown in Figure 8b, and 
clearly the singularities with smaller amplitudes are attracted to the points £ = ±1. 
When they get close enough to either of these points, an unperturbed finger continues 
to grow while the dimples remain near its base. If additional singularities are present 
that start much further away from the unit disk, they will arrive close to the unit disk 
at later times, producing new dimples near the tip of the finger. The endless presence 
of singularities streaming in from infinity can generate fingers with side-branching 
patterns seen frequently in experiments. 

In Figure 7a, the small amplitude of the singularity at £ 3 causes the formation 
of a small dimple on the side of a well-developed finger. In contrast, when Re E 3 
is the same order as Re E 1 and Re E 2 , the close approach of (3 to | £ |= 1 leads 
to a ‘tip-splitting’ event and the formation of two fingers which eventually compete. 
An example is given in Figure 9a for initial data of the form (23) with J\ — 2, J — 
3 and aj = -4/3 for j = 1,...,3. The initial singularity positions are chosen as 
(^(0) = — C 2 ( 0 ) = (2.5, 0.0), C3<0) = (0.0, 4.0) and the initial amplitudes as Ei((,\0) = 
(-1.0, 0.0), E 2 ( CO) = (-0.8, 0.0), £ 3 (C,0) = (-0.3, -0.16). We use N = 512 and 
At = 0.0005. As shown in Figure 9b, the singularities initially move radially towards 
the unit circle. As the singularity at £3 approaches \(\ — 1 , a dimple forms near the 
finger tip. As time advances, (3 gets closer to the interface and the dimple elongates 
into a large indentation, giving the appearance of two fingers. Eventually, the motion 
of the singularity is predominantly tangent to the circle, and the tangential velocity 
of £3 is such that the indentation is fixed in the laboratory frame. The direction the 
singularity moves around the circle determines which finger will dominate. In this case 
the singularity at (, 3 moves towards (,* = — 1. Since the point £ = — 1 corresponds to 
the bottom end of the interface, this motion has the effect of stretching the interface 
so that the indentation separating the fingers lies closer to the bottom end. The only 
way this can be done while also keeping the indentation fixed in the laboratory frame 
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is for the upper finger to grow significantly more than the lower one, as indicated in 
Figure 9a. Unfortunately, the simulation cannot be not run long enough to produce 
a clear outcome in the finger competition, because of the presence of the singularity 
in E J _ and G- at I/C3G). Nevertheless, we are able to compute long enough to make 
the trend in singularity motion clear. 

We conclude this section with a brief examination of some scenarios which may lead 
to a self-intersection of the interface. A self-intersection event is often referred to as a 
topological singularity; the possible formation of this type of singularity in Hele-Shaw 
flow and in other free surface flows is a topic of much current interest. A topological 
singularity occurs when the conformal map z((,t.) ceases to be univalent, i. e., when 
two points on the ( semi circle map to a single point. When this happens, either the 
more viscous or less viscous fluid region is divided into two disjoint sections. Bertozzi, 
Brenner, Dupont, & Kadanoff [7] and Goldstein, Pesci, & Shelley [12] have investigated 
possible topological singularity formation in Hele-Shaw flow with surface tension. They 
consider a particular geometry, consisting of a vertical cell with a thin layer of fluid 
resting on the bottom, chosen so that a variant of the lubrication approximation can 
be applied. Using this approximation, they concluded that in certain cases the top 
interface of the layer touches the bottom of the cell in finite time. However, their 
geometry and pressure conditions are significantly different from ours and it is unclear 
if their results can be extrapolated to our geometry, where there is a constant pressure 
gradient far ahead of the finger. 

In some situations, a loss of univalence is possible even when the singularities and 
zeroes of Z{((, t) in | ( |> 1 remain a finite distance from | ( |= 1. We used our code to 
search for such an event in the zero surface tension problem. Unlike commonly used 
boundary integral methods which run into resolution difficulties when the interface 
is about to pinch (see Baker & Shelley [3]), our numerical approach based on the 
conformal mapping function will not incur difficulties as long as the pinching is not 
accompanied by a singularity or zero of impinging on \ ( \= 1. 

Unfortunately, we were unable to find any occurance of this type of singularity 
in the situation where a fluid of zero viscosity displaces a viscous fluid. Of course, 
the formation of a topological singularity in this case cannot be completely ruled out 
from our limited examination and further study is required. When we reversed the 
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pressure gradient at infinity, however, so that the more viscous fluid on the right of 
Figure 1 displaces the fluid of negligible viscosity on the left, topological singularities 
were sometimes observed. It is well known from the Saffman & Taylor [36] analysis 
that a planar interface in this situation is stable and that any small deformation will 
reduce with time. This can be expected to be true even for most finite amplitude 
disturbances. However, our findings show that if the interface is highly deformed 
initially, pinching can occur. 

One such example is presented in Figure 10. This figure shows interfacial profiles 
before and after topological singularity formation for an initial value problem with 
two branch point singularities initially close to the unit disk. We picked data in 
of the form 0) = -2(1 - C 2 /C 2 )° r /( 7r 0 where a = -1.9 and ( s = (1.105,0.0); 
with this choice, zeros in z ^ are placed initially at infinityJI The change to liquid 
pushing air reverses the direction of the characteristic velocities (given by < 7 i),and 
the singularities move outward from the unit circle. Consequently, the evolution of 
the interface can be obtained without explicitly representing the singularities. We 
therefore set 2 - (2/7r)ln( = G+ + Go and solve (39) on the circle \(\ = e xv \ there is 
no need to evolve the radius of this circle. In the run we set N — 128 and At = 0.001. 

Our example shows that a loss of uni valency can occur in zero surface tension 
flows without concurrent singularities in 2 ^, and illustrates the ability of our numer- 
ical method to contend with such self intersections. The occurance of these singu- 
larities appears to be quite sensitive to initial conditions, with fatter initial fingers 
typically evolving to a flat sheet without pinching. The form of the solution after a 
self-intersection remains an open question. 

6 Discussion and Conclusions 

We have described a numerical method designed to track singularities present in the 
conformal map from the unit semi circle to the physical domain of Hele-Shaw flow in a 
channel. The method is restricted to those conformal maps that contain singularities 
of the form (23), and their initial location (^(0) and power a 3 must be given if they lie 

"Thus, the presence of surface tension will not significantly affect the evolution of the interface for 
the times shown. 
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inside the circle of radius 5(0). Furthermore, their amplitude 0) and G*(C, 0) must 
be specified, or at least their inner and outer components must be given on \(\ = 5(0) 
and |^*| = 1 respectively. This detailed information of the initial properties of the 
conformal map correspond to high precision knowledge of the initial interface location. 
Our numerical method then advances the conformal map, and hence the interface, 
numerically in a stable way. In other words, round-off errors do not contaminate the 
high precision specification of the interface location. 

Besides the restrictions on the initial conditions, there are two other limitations on 
our method. Singularities in inner analytic components of E 3 ((,t) and G(C, /) occur 
inside the unit disk at 1 /Cj(0* As Cj(G approaches the unit disk, these singularities also 
approach the boundary of the unit disk, and cause slow decay of the Fourier modes for 
E 3 ((,t) and G'(C\0* Consequently, we need many Fourier modes to ensure reasonable 
accuracy. Also, at present the method is limited to B = 0. When 5/0, zeros in 
are transformed into patterns of -4/3 singularities. Some of these can move toward 
the unit circle and afTect the shape of the interface at later times. Nevertheless, as long 
as Oj < -4/3 and none of the singularities formed out of initial zeros approach the 
unit disk closely, our results will be the correct limiting behavior as 5 — * 0. Perhaps 
more importantly, our method enables accurate 5 = 0 computations to be obtained 
for quite a general distributions of initial singularities, so that comparisons with 5 > 0 
solutions can be made. These kind of comparisons complement the asymptotic theory, 
and facilitate an understanding of the influence of small capillary effects. 

Despite the limitations, we find singularities induce interfacial structure that is 
typical of experiment observations when 5 is very small. In particular, two singulari- 
ties, placed on the real axis on either side of the unit circle, induce formation of a long 
finger. Singularities off the real axis induce small indentations on this finger if their 
amplitudes are small, giving the appearance of side-branches, or large indentations if 
their amplitudes are comparable to the ones on the real axis, giving the appearance 
of tip-splitting and finger competition. In general, we expect that a continual inward 
stream of singularities of all amplitudes can account for multiple branching and com- 
petition as observed experimentally. Although some aspects of interfacial evolution 
due to multiple singularities have been examined previously by Howison [19] using a 
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class of exact solutions for simple pole singularities, our algorithm allows a broader 
study to be undertaken for collections of isolated singularities of more general form. 

We plan to continue studies of the properties of the singularities and the interfa- 
cial structure their induce. For example, we wish to explore what role the complex 
amplitudes have on the trajectories of the singularities. We hope to find ways to rep- 
resent the singularities in better forms that may remove some of the limitations in our 
method, and we hope to find ways to capture the transformation of zeros when B is 
very small, but non-zero. The authors would like to acknowledge support from NASA 
grant NAG 3-1415 (G.B.), Department of Energy contract DE-FG02-92ER14270 (M. 
S. and S. T.), and an NSF Postdoctoral Fellowship (M. S.). S. T. was partially sup- 
ported by NASA grant NAS 1-18605 while in residence at The Institute for Computer 
Applications in Science and Engineering (ICASE). 
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A 




1. The unit semi circle in the ( plane is mapped into the viscous fluid region of the 
channel, with the circular arc being mapped to the interface. The points A,H, 
and C in the ( plane are mapped to the corresponding points in the channel. 


35 




36 


. The A = 1/2 SafFman-Taylor finger. In this and each ol the suMequcni ngures 
the viscous fluid region lies to the right of the curve. The profiles show the 
position of the interface from t = 0 to t = 3.0 in increments of .2. 


to 
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. A finger produced by two branch point singularities of strength rr = -4/5, 
located at symmetric positions on the real line. The profiles show the position 
of the interface from t = 0 to t = 0.8 in increments of .1, and at t = 0.85. In the 
last profile (t = 0.85) the singualrities are located at = ~( 2 = 1.0081. 
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. A finger produced by two branch point singularities of strength a = “4/3, 
located symmetrically on the real line. The profiles show the position of the 
interface from t = 0 to t = 1.1 in increments of .1. In the last profile ( t= 1 . 1 ) the 
singularities are located at £t = —£2 = 10054. 




39 


. Plot of log|E‘(A:)| versus -k for the computation in Figure 4. The plots range 
from t = 0 to t = 1.1 in increments of .1. Only odd values of k are plotted, 
since the values of E'(k) for even k are zero by symmetry (in the computed 
spectrum these values remain at round-ofT levels). The values of E l (k) for k > 1 
are identically zero. 
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a) Comparison of B = 0 and D = .0025 solutions for data witli two branch 
point singularities of strength o = -3.2 located on the real line, and with zeros 
of z ^ initially at infinity. The two solutions are indistinguishable at plotting 
resolution. The profiles are shown for < = 0 to / = .5 in increments of .1. 




41 


1025 (dashed line) solutions for 
shown for / = 0 to / = 1.1 in 
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b) The singularity trajectories corresponding to G(a). After reaching the unit 
circle, the (upper) complex singularity £3 moves in a counterclockwise direction. 
In the last profile (I = . 22 ), the singularities are located at = 1.1522, (2 = 
-1.1542, and <,3 = (-.9381, .3546). 
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pearence of side-branching. Here in 
t>= 0 to t = 1.2 in increments of .2. 



Re zeta 



(b) The singularities trajectories corresponding to 7(a). 
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a) An illustration of tip splitting followed by finger formation and competition 
due to a = -4/3 branch point singularities. The plots range from from t = 0 t< 
t = 0.6 in increments of .2, and from t = 0.6 to t = 0.7 in increments of .02. 
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b) The singularity trajectori 
motion of the singularities is 
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